library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(coda)
Rcpp::sourceCpp("/home/fra/Documents/GitHub/CommonAtomModel/Functions/CAM/newCAM.cpp")
source("/home/fra/Documents/GitHub/CommonAtomModel/Functions/CAM/newCAM.R")
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:purrr':
##
## compact
##
## Attaching package: 'reshape'
## The following objects are masked from 'package:plyr':
##
## rename, round_any
## The following object is masked from 'package:dplyr':
##
## rename
## The following objects are masked from 'package:tidyr':
##
## expand, smiths
## Registered S3 method overwritten by 'pryr':
## method from
## print.bytes Rcpp
##
## Attaching package: 'pryr'
## The following objects are masked from 'package:purrr':
##
## compose, partial
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2021 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
6 distribution osservate, 4 gruppi
u <- 30
x0 <- replicate(1, c(rnorm(u,3),rnorm( 3*u,10,.5)))
x1 <- replicate(1, c(rnorm(2*u,3),rnorm( 2*u,10,.5)))
x2 <- replicate(3, c(rnorm(2*u,-5),rnorm(2*u,10)))
x3 <- c(rnorm(2*u,3),rnorm( 2*u,10,.5),rnorm( 2*u,-5,.5),rnorm( 2*u,0,.5))
yd <- (c(x0,x1,x2,x3))
yg <- rep(1:6,c(rep(4*u,5),8*u))
plot(yd,col=yg)
nb <- 10000
m1 <- CAM(y_obser = yd,
y_group = yg,
K0 = 20,
L0 = 20,
prior = list(
# hyperparameters NIG
m0= mean(yd),
k0=1/var(yd),
a0=3, b0=1,
# hyperparameters alpha and beta
a_alpha=3, b_alpha =3,
a_beta =3, b_beta = 3),
nsim = nb,
burn_in = nb,
thinning = 1,
verbose = F,
fixedAB = T,
kappa=0.5,
cheap=F,
seed=12345,
sigma.prop.beta = 1)
m2 <- CAM(y_obser = yd,
y_group = yg,
K0 = 20,
L0 = 20,
prior = list(
# hyperparameters NIG
m0= mean(yd),
k0=1/var(yd),
a0=3, b0=1,
# hyperparameters alpha and beta
a_alpha=3, b_alpha =3,
a_beta =3, b_beta = 3),
nsim = nb,
burn_in = nb,
thinning = 1,
verbose = F,
fixedAB = F,
kappa=0.5,
cheap=F,
seed=123456,sigma.prop.beta = 1,
conditional.beta = F
)
m3 <- CAM(y_obser = yd,
y_group = yg,
K0 = 20,
L0 = 20,
prior = list(
# hyperparameters NIG
m0= mean(yd),
k0=1/var(yd),
a0=3, b0=1,
# hyperparameters alpha and beta
a_alpha=3, b_alpha =3,
a_beta =3, b_beta = 3),
nsim = nb,
burn_in = nb,
thinning = 1,
verbose = F,
fixedAB = F,
kappa=0.5,
cheap=F,
seed=123456,sigma.prop.beta = 1,
conditional.beta = T
)
plot(as.mcmc(m1$A_DP),main="M1")
plot(as.mcmc(m2$A_DP),main="M2")
plot(as.mcmc(m3$A_DP),main="M3")
par(mfrow=c(1,2))
acf(m2$A_DP)
pacf(m3$A_DP)
par(mfrow=c(3,1))
plot(as.mcmc(m1$B_DP),main="M1")
plot(as.mcmc(m2$B_DP),main="M2")
plot(as.mcmc(m3$B_DP),main="M3")
par(mfrow=c(1,2))
acf( m2$B_DP)
pacf(m3$B_DP)
par(mfrow=c(1,3))
matplot(m1$NN.cz,type="l")
matplot(m2$NN.cz,type="l")
matplot(m3$NN.cz,type="l")
par(mfrow=c(1,1))
plot(density(m1$NN.cz[,1],type="l") ,col=1,lwd=2)
## Warning: In density.default(m1$NN.cz[, 1], type = "l") :
## extra argument 'type' will be disregarded
lines(density(m2$NN.cz[,1],type="l"),col=2,lwd=2)
## Warning: In density.default(m2$NN.cz[, 1], type = "l") :
## extra argument 'type' will be disregarded
lines(density(m3$NN.cz[,1],type="l"),col=3,lwd=2)
## Warning: In density.default(m3$NN.cz[, 1], type = "l") :
## extra argument 'type' will be disregarded
plot(density(m1$NN.cz[,2],type="l") ,col=1,lwd=2)
## Warning: In density.default(m1$NN.cz[, 2], type = "l") :
## extra argument 'type' will be disregarded
lines(density(m2$NN.cz[,2],type="l"),col=2,lwd=2)
## Warning: In density.default(m2$NN.cz[, 2], type = "l") :
## extra argument 'type' will be disregarded
lines(density(m3$NN.cz[,2],type="l"),col=3,lwd=2)
## Warning: In density.default(m3$NN.cz[, 2], type = "l") :
## extra argument 'type' will be disregarded
r1=PSM(m1$Z_j)
image(r1)
a1=mcclust.ext::minVI(r1,method = "greedy")
plot(yd,col=rep(a1$cl,table(yg)))
r2=PSM(m2$Z_j)
image(r2)
a2=mcclust.ext::minVI(r2,method = "greedy")
plot(yd,col=rep(a2$cl,table(yg)))
r3=PSM(m3$Z_j)
image(r3)
a3=mcclust.ext::minVI(r3,method = "greedy")
plot(yd,col=rep(a3$cl,table(yg)))
library(reshape2)
##
## Attaching package: 'reshape2'
## The following objects are masked from 'package:reshape':
##
## colsplit, melt, recast
## The following object is masked from 'package:tidyr':
##
## smiths
xx <- seq(-10,15,by = .1)
LLL <- list()
LLL1 <- posterior.densities(m1,howfarback = 100,normalize = T,xx = seq(-10,15,by = .10))
## 123456
D <- melt(LLL1)
p1 <- ggplot(D)+
geom_line(aes(x=Var2,y=value,group=interaction(Var1,L1),col=factor(L1)),alpha=.25)+
theme_bw()+scale_color_viridis_d()
LLL2 <- posterior.densities(m2,howfarback = 100,normalize = T,xx = seq(-10,15,by = .10))
## 123456
D <- melt(LLL2)
p2 <- ggplot(D)+
geom_line(aes(x=Var2,y=value,group=interaction(Var1,L1),col=factor(L1)),alpha=.25)+
theme_bw()+scale_color_viridis_d()
LLL3 <- posterior.densities(m3,howfarback = 100,normalize = T,xx = seq(-10,15,by = .10))
## 123456
D <- melt(LLL3)
p3 <- ggplot(D)+
geom_line(aes(x=Var2,y=value,group=interaction(Var1,L1),col=factor(L1)),alpha=.25)+
theme_bw()+scale_color_viridis_d()
library(patchwork)
##
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
##
## area
p1/p2/p3
par(mfrow=c(3,1))
r1 <- apply(m1$post.dens,c(1,2),mean)
matplot(r1,type="l",main="M1")
r2 <- apply(m2$post.dens,c(1,2),mean)
matplot(r2,type="l",main="M2")
r3 <- apply(m3$post.dens,c(1,2),mean)
matplot(r3,type="l",main="M3")